function sb = solve_backward(param, fp,fp_parfor)

alpha = param(1:9);
TFP = param([10:11 end]);
delta = param([12:14 (end-1)]);
gamma = param(15:19);

delta(2) = delta(2) + param(20)*(fp.mean_income_by_neighborhoods(fp_parfor.j) - fp.mean_income_by_neighborhoods(4) );

beta = 0.95;

estim_policy0=cell(fp.periods-1,2);
estim_policy1=cell(fp.periods-1,2);

estim_value  =cell(fp.periods-1,2);
estim_value0 =cell(fp.periods-1,2);
estim_value1 =cell(fp.periods-1,2);

%Grid points;
grid_peers=fp_parfor.grid_peers;
inv0 = fp_parfor.inv0;
inv1 = fp_parfor.inv1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Create Grid;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

hk = fp_parfor.grid ;

if fp_parfor.i==1
    %Initial guess in the solution algorithm;    
    Hc_tp1_guess{1} = max(fp_parfor.distrib_skills(:,2)).*ones( length(inv1{1}), length(hk{1}) , length(grid_peers{1}) , 2 );
    Hc_tp1_guess{2} = max(fp_parfor.distrib_skills(:,3)).*ones( length(inv1{2}), length(hk{2}) , length(grid_peers{2}) , 2  );

else
    % Previous distribution in the fixed-point algorithm;    
    Hc_tp1_guess = fp_parfor.Hc_tp1_guess_prime ;

end


for t=fp.periods-1:-1:1

    grid_p=grid_peers{t};
    grid_p = sort(grid_p);


    if t==fp.periods-1     %Solve last period

        grid_inv0 = inv0{t};
        grid_inv1 = inv1{t};

        [I0, hc0, H0, College0] =ndgrid(grid_inv0, hk{t}, grid_p, [0,1] );
        [I1, hc1, H1, College1] =ndgrid(grid_inv1, hk{t}, grid_p, [0,1] );

        log_Skills_tp1_0 = log(technology(t,alpha, hc0, I0, 0, H0, College0, TFP) ) ;
        log_Skills_tp1_1 = log(technology(t,alpha, hc1, I1, 1, H1, College1, TFP) ) ;

        U0 = delta(1).*log(hc0).*(1 - delta(3).*0) + 1.*log(1-I0)  -delta(2).*0 + beta.*delta(1).*log_Skills_tp1_0;
        U1 = delta(1).*log(hc1).*(1 - delta(3).*1) + 1.*log(1-I1)  -delta(2).*1 - delta(4).*(College1==1) + beta.*delta(1).*log_Skills_tp1_1;
        
        %Find the optimal parental investments by parenting style;
        [Value0,ind0] = max(U0,[],1);
        [Value1,ind1] = max(U1,[],1);


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Computing Value e Policy Functions;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;


        for c = 1:1:2 %loop over education (c=1 non-college, c=2 college)

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
            %Value Function;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
            Value0_actual(:,:,:,c,t) = Value0(:,:,:,c);
            Value1_actual(:,:,:,c,t) = Value1(:,:,:,c);


            %Collect State Variables;
            hc_s = hc0(1,:,:, c);
            H_s  = H1(1,:,:, c);
            State_var = [ hc_s(:), H_s(:) ];

            %Value for P=0;
            value0 = Value0(:,:,:,c);           
            Y0 = value0(:);
            estim_value0{t, c} = compute_value_function(Y0, State_var);    

            %Value for P=1;
            value1 = Value1(:,:,:,c);
            Y1 = value1(:);
            estim_value1{t, c} = compute_value_function(Y1, State_var);    

             %Expected Value over preference shocks;
            Y = log( exp(Y0) + exp(Y1));
            estim_value{t,c}  = compute_value_function(Y, State_var);            

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Policy Function (Investments)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            %Store the optimal choice for each state (by P);
            Inv0 = I0(ind0);
            Inv1 = I1(ind1);

            inv0_temp =  Inv0(:,:,:,c);            
            Y0 = inv0_temp(:);
            estim_policy0{t,c} = compute_policy_function(Y0, State_var);

            inv1_temp =  Inv1(:,:,:,c);                        
            Y1 = inv1_temp(:);            
            estim_policy1{t,c} = compute_policy_function(Y1, State_var);


        end

    else %Before Last period;

        grid_inv0 = inv0{t};
        grid_inv1 = inv1{t};


        [I0, hc0, H0, College0] =ndgrid(grid_inv0, hk{t}, grid_p, [0,1] );
        [I1, hc1, H1, College1] =ndgrid(grid_inv1, hk{t}, grid_p, [0,1] );

        %Distribution coming from the Solve_Forward;
        set_kids2 = fp_parfor.distrib_skills(:,t+1);
        set_PS2   = fp_parfor.distrib_PS(:,t);

        %Iteration
        set_kids = Hc_tp1_guess{t};
        set_kids1 = set_kids(:);

        %Random graph generation;
        [H1_tp1, H2_tp1] = ndgrid(set_kids1, set_kids2);

        set_PS_0 = zeros(length(set_kids1),1);
        set_PS_1 = ones(length(set_kids1),1);

        [PS1_0, PS2_0] = ndgrid(set_PS_0, set_PS2);
        [PS1_1, PS2_1] = ndgrid(set_PS_1, set_PS2);

        m_p0 = mean_peers(H1_tp1,H2_tp1,PS1_0,PS2_0,gamma,fp,set_kids2,fp_parfor,t);
        H_tp1_0 = m_p0.H_prime ;

        m_p1 = mean_peers(H1_tp1,H2_tp1,PS1_1,PS2_1,gamma,fp,set_kids2,fp_parfor,t);
        H_tp1_1 = m_p1.H_prime ;

        %Evolution of peer qualily by parenting style;
        H_prime_0 = reshape(nanmean(log(H_tp1_0),2),[size(hc0,1),size(hc0,2),size(hc0,3),size(hc0,4)]);
        H_prime_1 = reshape(nanmean(log(H_tp1_1),2),[size(hc1,1),size(hc1,2),size(hc1,3),size(hc1,4)]);

        %Evolution of a child's skills by parenting style;        
        log_Skills_tp1_0 = log(technology(t,alpha, hc0, I0, 0, H0, College0, TFP) ) ;
        log_Skills_tp1_1 = log(technology(t,alpha, hc1, I1, 1, H1, College1, TFP) ) ;

        %Calculate Altruism;

        %Calculate Altruism based on Expected Value of Child's Utility (P=0);
        temp_utility_child0 = nanmean( m_p0.U_links1,3);
        utility_child0 = nansum( temp_utility_child0,2);
        Altruism_0 = reshape( utility_child0 ,[size(hc0,1),size(hc0,2),size(hc0,3),size(hc0,4)]);

        %Calculate Altruism based on Expected Value of Child's Utility (P=1);        
        temp_utility_child1 = nanmean( m_p1.U_links1,3);
        utility_child1 = nansum( temp_utility_child1,2);
        Altruism_1 = reshape( utility_child1 ,[size(hc1,1),size(hc1,2),size(hc1,3),size(hc0,4)]);



        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Solving the parent's problem;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        for c = 1:1:2 %loop over education (c=1 non-college, c=2 college)

            %Continuation value by parenting style (given parent education);
            Cont_value0(:,:,:,c) = Interp_function(estim_value{t+1,c},log_Skills_tp1_0(:,:,:,c),H_prime_0(:,:,:,c));
            Cont_value1(:,:,:,c) = Interp_function(estim_value{t+1,c},log_Skills_tp1_1(:,:,:,c),H_prime_1(:,:,:,c));

            %Bellman Equation by Parenting Style and Education of parents (c);
            U0(:,:,:,c) = fp.lambda.*delta(1).*log(hc0(:,:,:,c)).*(1 - delta(3).*0) + (1-fp.lambda).*Altruism_0(:,:,:,c) + 1.*log(1-I0(:,:,:,c))  -delta(2).*0 + beta.*Cont_value0(:,:,:,c);
            U1(:,:,:,c) = fp.lambda.*delta(1).*log(hc1(:,:,:,c)).*(1 - delta(3).*1) + (1-fp.lambda).*Altruism_1(:,:,:,c) + 1.*log(1-I1(:,:,:,c))  -delta(2).*1 - delta(4).*(c==2) + beta.*Cont_value1(:,:,:,c);            

        end

        %Find the optimal parental investments by parenting style;
        [Value0,ind0] = max(U0,[],1);
        [Value1,ind1] = max(U1,[],1);

        temp_value0 = repmat( Value0 , [size(I0,1), 1, 1 , 1] );
        temp_value1 = repmat( Value1 , [size(I1,1), 1, 1 , 1] );


        %Optimal behavior about parenting style;
        I_guess0 = repmat(I0(ind0), [size(I0,1), 1, 1 , 1] );
        I_guess1 = repmat(I1(ind1), [size(I1,1), 1, 1 , 1] );

        Prob_PS = (1./(1+1./exp(temp_value1 -temp_value0)));

        %Updating the guess for next-period peers skills according to optimal behavior of parents;
        Hc_tp1_guess{t} = (1-Prob_PS).*technology(t,alpha, hc1, I_guess0, 0, H1, College0, TFP) + Prob_PS.*technology(t,alpha, hc1, I_guess1, 1, H1, College1, TFP);


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Computing Value e Policy Functions;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        for c = 1:1:2 %loop over education (c=1 non-college, c=2 college)


            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
            %Value Function;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
            Value0_actual(:,:,:,c,t) = Value0(:,:,:,c);
            Value1_actual(:,:,:,c,t) = Value1(:,:,:,c);


            %Collect State Variables;
            hc_s = hc0(1,:,:, c);
            H_s  = H1(1,:,:, c);
            State_var = [ hc_s(:), H_s(:) ];

            %Value for P=0;
            value0 = Value0(:,:,:,c);           
            Y0 = value0(:);
            estim_value0{t, c} = compute_value_function(Y0, State_var);    

            %Value for P=1;
            value1 = Value1(:,:,:,c);
            Y1 = value1(:);
            estim_value1{t, c} = compute_value_function(Y1, State_var);    

             %Expected Value over preference shocks;
            Y = log( exp(Y0) + exp(Y1));
            estim_value{t,c}  = compute_value_function(Y, State_var);


            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Policy Function (Investments)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            %Store the optimal choice for each state (by P);
            Inv0 = I0(ind0);
            Inv1 = I1(ind1);

            inv0_temp =  Inv0(:,:,:,c);            
            Y0 = inv0_temp(:);
            estim_policy0{t,c} = compute_policy_function(Y0, State_var);

            inv1_temp =  Inv1(:,:,:,c);                        
            Y1 = inv1_temp(:);            
            estim_policy1{t,c} = compute_policy_function(Y1, State_var);

        end %end loop for parental education "c"

    end  %End "if last period or not";

end %End loop for "t";

sb.Hc_tp1_guess_prime = Hc_tp1_guess;

sb.estim_policy0 = estim_policy0;
sb.estim_policy1 = estim_policy1;
sb.estim_value0  = estim_value0;
sb.estim_value1  = estim_value1;
sb.Value0_actual = Value0_actual;
sb.Value1_actual = Value1_actual;
